#!/usr/bin/perl
# 2013/05/29 Calcualte 3kb(average binding distance) from summit of peak
# 2013/05/08 Categorize association types as binding site associated or non associated

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -r [peak sumit file] -o [output filtered map file]\n";
}

my %opt;
getopts("i:o:r:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};
my $FILE_peaks = $opt{r};



#---------------------------------------
# read peaks file
#---------------------------------------
my %peaks;
my $fh_peaks = IO::File->new($FILE_peaks) or die "cannot open $FILE_peaks: $!";
while($_ = $fh_peaks->getline()){
	s/\r?\n//;
	my ($chr, $start, $end, $name, $pval) = split /\t/;

	# take surrounded 3kb
	$start -= 1500;
	$end = $start + 2999;

	push @{$peaks{$chr}}, [$start, $end, $name];
}
$fh_peaks->close();

#---------------------------------------
# connect binding information
#---------------------------------------
my $COUNTER_bind_bind = 0;
my $COUNTER_bind_other = 0;
my $COUNTER_other_other = 0;

my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";

# id chr1 loc1 direction1 quality1 hind1 hind_loc1 chr2 loc2 direction2 quality2 hind2 hind_loc2 type(unique or repeat) number(alignable locations) locations
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	### check reads is overlapped with peaks or not
	my ($peakName1, $peakLoc1) = &checkPeaks($chr1, $loc1);
	my ($peakName2, $peakLoc2) = &checkPeaks($chr2, $loc2);

	### count
	my $combination = 0;
	if($peakName1 ne "NA"){
		$combination++;
	}
	if($peakName2 ne "NA"){
		$combination++;
	}

	if($combination == 0){
		$COUNTER_other_other++;
	}elsif($combination == 1){
		$COUNTER_bind_other++;
	}else{
		$COUNTER_bind_bind++;
		$fh_out->print(join("\t", $id, $chr1, $loc1, $direction1, $mapQ1, $peakName1, $peakLoc1, $chr2, $loc2, $direction2, $mapQ2, $peakName2, $peakLoc2, $type, $num, $locations) . "\n");
	}
}
$fh_map->close();
$fh_out->close();

my $SUM = $COUNTER_bind_bind + $COUNTER_other_other + $COUNTER_bind_other;
printf "Binding - Binding: %d (%.3f)\n", $COUNTER_bind_bind, $COUNTER_bind_bind / $SUM * 100;
printf "Binding - Other : %d (%.3f)\n", $COUNTER_bind_other, $COUNTER_bind_other / $SUM * 100;
printf "Other - Other : %d (%.3f)\n", $COUNTER_other_other, $COUNTER_other_other / $SUM * 100;

sub checkPeaks{
	my ($chr, $pos) = @_;
	foreach my $m(@{$peaks{$chr}}){
		my ($peakS, $peakE, $name) = @$m;
		my $location = $peakS . '-' . $peakE;
		if($peakS <= $pos and $pos <= $peakE){
			return ($name, $location);
		}
	}
	return ("NA", "NA");
}

